function [rho,rhou,E] = LaxFriedrichs3(n,v)
%用LaxFridrichs格式计算Euler方程的初边值问题3
%参数说明：
%   输入:
%   T：输出时刻
%   n：单位长度上的网格数量。
%   v:\tau / h
%   varargin：输出选项。
%   输出：
%   rho,rhou,E:相应变量在区T时刻的值.
    T = 0.0381;
    T3=[0.010,0.016, 0.026, 0.028, 0.030, 0.032, 0.034, 0.038,100];
    
    h = 1/n;
    tau = h*v;
    step = floor(T/tau);
    %初始化
    rho = zeros(n+1,2);
    rhou = zeros(n+1,2);
    E = zeros(n+1,2);

    %初始条件
    rho(:,1)=1;
    rhou(:,1)=0;
    a = floor(0.1*n);
    b = floor(0.9*n);

    E(1:a,1)=2.5*10^3;
    E(a+1:b,1)=2.5*10^-2;
    E(b+1:n+1,1)=2.5*10^2;
    
    %边界条件
    rho(1,2)=rho(1,1);
    rho(n+1,2)=rho(n+1,1);
    E(1,2)=E(1,1);
    E(n+1,2)=E(n+1,1);

    
    frho=zeros(1,n+1);
    frhou=zeros(1,n+1);
    fE=zeros(1,n+1);
    %迭代
    ti = 2;
    t = 1;
    printtime = 1;
    for s = 2:step+1
        x = t;
        t = ti;
        ti = x;

        %计算f_j
        rhou2 =(rhou(1:n+1,ti).^2./rho(1:n+1,ti));
        p = 0.4*(E(1:n+1,ti)-rhou2/2);
        frho(1:n+1)=rhou(1:n+1,ti);
        frhou(1:n+1)=rhou2+p;
        fE(1:n+1)=rhou(1:n+1,ti)./rho(1:n+1,ti).*(E(1:n+1,ti)+p);
        
        for i=2:n
            rho(i,t)=(rho(i-1,ti)+rho(i+1,ti))/2-v/2*(frho(i+1)-frho(i-1));
            rhou(i,t)=(rhou(i-1,ti)+rhou(i+1,ti))/2-v/2*(frhou(i+1)-frhou(i-1));
            E(i,t)=(E(i-1,ti)+E(i+1,ti))/2-v/2*(fE(i+1)-fE(i-1));
        end
        now = (s-1)/step*T;
        if now>T3(printtime)
            printtime = printtime+1;
            ax = 0:1/n:1;
            plot(ax,rho(:,t));
            pause(0.1)
            title(['Question 3 \rho at time ',num2str(now)]);
            h = gcf;
            hgexport(h,['report/LF/q3_time_',num2str(now),'.eps']);
        end
    end
    rho = rho(:,t);
    rhou = rhou(:,t);
    E = E(:,t);
end